A Multi-level Algorithm for Quantum-impurity Models 
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A continuous-time path integral Quantum Monte Carlo method using the directed-loop algorithm 
is developed to simulate the Anderson single-impurity model in the occupation number basis. Al- 
though the method suffers from a sign problem at low temperatures, the new algorithm has many 
advantages over conventional algorithms. For example, the model can be easily simulated in the 
Kondo limit without time discretization errors. Further, many observables including the impurity 
susceptibility and a variety of fermionic observables can be calculated efficiently. Finally the new 
approach allows us to explore a general technique, called the multi-level algorithm, to solve the sign 
problem. We find that the multi-level algorithm is able to generate an exponentially large number 
of configurations with an effort that grows as a polynomial in inverse temperature such that config- 
urations with a positive sign dominate over those with negative signs. Our algorithm can be easily 
generalized to other multi-impurity problems. 

PACS numbers: 02.70.Ss, 05.30.Fk, 72.15.Qm 



I. INTRODUCTION 

Developing efficient algorithms to solve problems 
in statistical mechanics involving strongly correlated 
fermions is an important area of research in computa- 
tional physics. Such problems are often afflicted with sign 
problems and are difficult to handle using Monte Carlo 
techniques. The conventional approach is to integrate the 
fermions out and write the problem in terms of a statis- 
tical mechanics of a bosonic system where the fermionic 
physics is hidden in the Boltzmann weight as a deter- 
minant of a matrix 0, 01 ■ In some interesting cases the 
determinant is positive definite which makes it possible 
to design Monte Carlo methods for solving the problem. 
This approach is often referred to as the determinantal 
Monte Carlo method. Unfortunately, the method fails 
in many cases since the fermion determinant can often 
be negative. Even when the determinant is guaranteed 
to be positive the algorithms can be inefficient in certain 
interesting range of parameters. Determinantal Monte 
Carlo methods can often be formulated only when the 
partition function is written as a path integral in dis- 
crete Euclidean time, which leads to time-discretization 
errors. Although such errors are controllable through ex- 
trapolation techniques they can be time consuming. All 
these difficulties make it important to explore alternative 
algorithms. 

Recently, another approach to fermionic physics in 
which fermionic partition functions are written in the 
occupation number basis has gained popularity 0, 
Although in this approach one encounters sign problems 
it is sometimes possible to design efficient algorithms in 
regions of parameter space where the sign problems are 
mild. In certain cases this approach also leads to novel 
solutions to the sign problems 0, which in turn lead to 
algorithms which are much more efficient than conven- 
tional ones. This has lead recently, for instance, to the 
first successful confirmation of the Kosterlitz-Thouless 
behavior in a fermionic model lal. In this article, we ex- 



plore a new algorithm to study the physics of electrons in 
a partially filled band interacting with a few impurities. 

Quantum impurity models are, of course, classics of 
condensed matter many-body physics 0- Interest in 
them has been reawakened in recent years because of 
developments from two completely different points of 
view. On the one hand, quantum dots in semiconductor 
heterostructures allow for the creation of tunable quan- 
tum impurities which can be studied individually and 
with exquisite precision 8]. Our particular interest is 
in studying the effects of mesoscopic fluctuations on the 
many-body physics On the other hand, the study 
of strongly correlated electron systems away from half 
filling - systems which continue to be investigated inten- 
sively - leads naturally to fermionic quantum impurity 
models through the dynamical mean theory approxima- 
tion We plan to use our method in this connection 
in the future. Since the Anderson single impurity model 
is the simplest in this class 11, J^], we focus on it; it is 
straight forward to extend the method to include more 
impurities. 

The Hamiltonian of the Anderson impurity model 
which we consider is 

k,o- <T 

+ 5I^k(cLrf. + 4ck.) + C/4dt4dx, (1) 
k,o- 

where the first term represents the free band electron 
energy levels, the second term is the impurity energy, 
the third represents the hopping amplitudes between free 
electron states and the impurity, and the last term is the 
repulsive coulomb interaction term at the impurity site. 
We assume that —D<e\^<D where 2D \s the band- 
width. This model was introduced over forty years ago by 
Anderson Tl'l to study the affects of impurity spins em- 
bedded in metals. Today, this model plays a very impor- 
tant role in understanding a variety of condensed matter 
systems . The problem can be solved analytically for 
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the case of a constant hopping ampHtude and a constant 
density of energy levels with an infinite bandwidth, us- 
ing the Bethe ansatz [T^ IT3 | . In the limit of large U one 
can relate this model to the famous Kondo model [T^ . 
After his discovery of renormalization, Wilson used this 
model to illustrate the numerical renormalization group 
program which is a powerful method to solve this 
problem. It is now well understood that the low temper- 
ature properties of this model require a non-perturbative 
approach. 

More than a decade ago Hirsch and Fye developed a de- 
terminantal Monte-Carlo algorithm to study this model 
0- Such an approach is necessary to solve the Hamilto- 
nian in Eq.QJ with general parameters. In this method 
the partition function for a given temperature is rewrit- 
ten as a statistical mechanics of a model on a discrete 
temporal lattice. An auxiliary Ising variable for the im- 
purity is introduced as a function of time in order to con- 
vert the Hamiltonian into a problem in which fermions 
are free but interact with the Ising spin. It is then pos- 
sible to integrate the fermions completely to write the 
problem as a statistical mechanics of the Ising spin. The 
Boltzmann weight can be written in terms of certain 
Green functions which can be computed easily. Unfortu- 
nately, the algorithm slows down as the temporal lattice 
spacing is taken to zero for a fixed temperature. It is 
also difficult to approach the large U limit since the dis- 
cretization error becomes significant. Finally, some ob- 
servables are difficult to evaluate. One famous example 
is the impurity su scep tibility which is known to contain 
large fluctuations |21|. 

Here we explore an alternate approach, by writing the 
partition function in the occupation number basis in con- 
tinuous Euclidean time. Since fermionic occupations con- 
sist of two states one can use the recently developed 
dircctcd-loop algorithm for quantum spin systems |l6| 
in continuous time Unlike the Hirsch and Fye al- 

gorithm, in our method one can only deal with a finite 
number of energy levels, but the discretization error in 
the Euclidean time direction can be eliminated. This al- 
lows us to simulate a large value of U with little effort. 
In the occupation number basis we can also easily calcu- 
late observables such as the average occupation numbers, 
and the local susceptibility. Moreover, as we will discuss 
later we are able to calculate the impurity susceptibilities 
efficiently. 

Inspite of these advantages, the Boltzmann weights of 
our configurations can be negative due to the fermion 
permutation sign. Thus, our method suffers from a sign 
problem. However, we find that the sign problem is 
rather mild down to the Kondo temperature. Inter- 
estingly, our approach also allows us to explore a new 
technique, called the multi-level algorithm, to solve the 
sign problem. This technique was recently used in lattice 
QCD in determi ning the string tension between quarks 
and anti-quarks [l9l |. A similar technique was also ex- 
plored in ha. Since the multi-level technique is a general 
method, our method allows us to study the usefulness of 



this approach to solve a class of fermion sign problems. 
Here we show that the multi-level technique is indeed 
useful in alleviating the sign problem. We were able to 
estimate signs of the order of 10~® using this approach. 

The paper is organized as follows: In section II we in- 
troduce the Monte Carlo algorithm for updating the An- 
derson impurity model in the occupation number basis. 
In section HI we explain the multi-level algorithm and 
show how one can use it to calculate the average signs 
efficiently. In section IV we discuss how one can measure 
observables in our method. We discuss the calculation of 
the impurity susceptibility and how our method allows us 
to compute it efficiently. Section V contains some results 
from the algorithm. Section VI contains our conclusions. 



II. THE DIRECTED LOOP ALGORITHM 

In this section we construct a Monte Carlo method 
to calculate quantities for the Anderson impurity model 
described by the Hamiltonian of Eq.(^. We begin by 
rewriting the partition function, Z — Tr e~^^ at tem- 
perature T = 1/(3, as a path integral in Euclidean time. 
This is accomplished by introducing M(= P/t) imagi- 
nary time slices and writing 
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C i=0 



where Ci represent the electron states on the i-th time 
slice in the occupation number basis. Since we are evalu- 
ating the trace we must have Co = Cm- The true parti- 
tion function in continuous time is obtained in the limit 
of large M and small r at fixed (3. We can define 



M-l 



W[C]ct[C] = n < C^\e-^"\C,+i >, 



(3) 



where W{C) is the magnitude and a[C] the sign of the 
Boltzmann weight. Then the partition function can be 
written as 



Z = J2w[C]a[C]. 



(4) 



In the Monte Carlo method each space-time configuration 
C is generated stochastically with probability 



PiC) = 



W[C] 



(5) 



Ignoring the sign, each configuration of fermion occupa- 
tion numbers is analogous to that of a configuration of a 
quantum spin-half particles. Hence we can use an exten- 
sion of the directed loop algorithm discussed in to 
update the space-time occupation number configurations 
C directly in continuous time. A simple way to construct 
such an update is to construct the update rules for finite 



FIG. 1: When the directed loop enters a vertex (shown on 
the left) one can produce many exit paths (some of them are 
shown on the right after the flip in the occupation numbers). 
Circles indicate electron levels in the band and squares rep- 
resent the impurity (filled symbols indicate occupied sites). 
Each figure only shows a few relevant electron sites. For the 
Anderson model discussed in the text (c) and (d) are forbid- 
den to order t. Both (a) and (b) are allowed but we can 
choose the weight of (a) to be zero while satisfying detailed 
balance. 

M and then take the hmit of infinite AI. Below we de- 
scribe our rules. 

A special feature of the Anderson model is that elec- 
tron hopping must include the impurity and any of the 
band electron sites. However, since all the band electron 
sites are involved, in our algorithm a "vertex" , in the lan- 
guage of |l6j| , is change in the configuration between two 
adjacent time slices, Ci and Ci+i. If Ci ~ C^+i then one 
has a diagonal configuration with weight 

< Cle^^-^ICj+i > = l-T^etrik^ 

k,cr 

a 

up to first order in r. Since we will take the continuum 
time limit this is sufficient. If C,; ^ C^+i then the two 
configurations can only differ in the occupation numbers 
of either the spin up or the spin down impurity fermion 
and the corresponding spin of one of the band fermion 
levels with momentum say k. Further both spins cannot 
hop simultaneously! Thus, if these constraints are met 
then 

<Q|e-^^|a+i >=t|14|, (6) 

otherwise the matrix element is zero and is disallowed. 
Our directed loop update begins by choosing a point, 
with 50% probability on the impurity and 50% probabil- 
ity on the other sites chosen randomly on one of the spin 
layers. Then with probability half the path enters the 
vertex either in the positive time direction or the neg- 
ative time direction. Using a set of rules that governs 
the exit of the update path at each "vertex" given the 
entrance of the path, the loop grows until it finds the 
starting point, where the update ends. The occupation 
numbers are changed while the loop is being constructed. 

Since the Hamiltonian commutes with the total spin 
operator and the total fermion number operator, the 
rules that generate the loop must satisfy the appropriate 
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FIG. 2: Assignments of directed-loop segments. All possible 
vertices in which the entrance site is occupied are shown. Ver- 
tex E was considered in Fig.0 However, unlike Fig. all the 
exit paths are shown together without the flip in occupation 
numbers. 

conservation laws. In addition, since we will ultimately 
take the r ^ limit, processes of order should not 
be considered. Thus, only one hopping can occur at a 
vertex. Using these constraints one can determine all the 
allowed processes in a given vertex. This is illustrated 
in Fig. n in which the directed loop enters the vertex 
through a filled electron state in the band with a definite 
spin. On the right, four occupation number conserving 
exit paths ((a),(b),(c),(d)) are shown. It is easy to see 
that, up to order r, only (a) and (b) are allowed. Once all 
the possible paths are determined, all the processes are 
given probability weights such that they satisfy detailed 
balance. For example in this case the bounce weight (a) 
can be chosen to be zero, so that the continuation process 
(b) occurs with probability one. For vertices where the 
path enters an occupied site, all possible loop segment 
assignments are shown in Fig. |21 The example consid- 
ered in Fig. n is shown as vertex E in Fig. [21 Let us now 
discuss the other vertices. 

First, consider a vertex where a hopping occurs at site 
k and the loop enters the vertex at site fc (vertex A in 
Fig. 121). There are a total of -I- 1 different possible 
exits in this vertex where N is the total number of band 
energy levels. Let A\^q be the weight for the path to exit 
at another band electron level q. We choose 

min[|14|,|y,|] 
= ^ N+1 ■ 

The probability for this process will then be given by 
Pkq = ^feg/(r| VkD- Since Akq is symmetric in k and q 
detailed balance is satisfied. The remaining probability 
must be the probability for the loop to exit at the im- 
purity. Since there are two possible paths, for simplicity 
the weight for each of these processes is chosen to be 

Akd = l{r\Vk\-J2A,k). (8) 

q^k 

Then, the probability for these process is Afe£;/T|Vk|. The 
factor A'' -|- 1 guarantees that all of Aq^. , Adk are posi- 
tive numbers. Our choice of weights is such that in the 
case where all |Vk| are the same all possible processes are 
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equally likely. Note that the loop does not bounce back 
at this vertex, and there is a unique direction for the loop 
at each of the band levels. 

In the case of vertex B the incoming path is on the 
impurity. The outgoing path can be at one of the elec- 
tron band levels with momentum k which has the weight 
Adk which will be chosen to be equal to Akd to satisfy 
detailed balance. In this case there is a possibility for 
the loop to continue forward or to bounce back. In order 
to fix the weights of these processes one compares the 
weight of the original vertex with that of the vertex ob- 
tained if the loop continues to go forward. If the forward 
continuation produces a vertex of the smaller weight, the 
bounce weight is chosen to be the absolute value of dif- 
ference of the two weights; otherwise the bounce weight 
is zero. There are two cases to be considered: for the 
case in which the impurity site contains an electron with 
the opposite spin (to the spin at which the path is being 
constructed), the bounce weight is \T{ed + U)\, and in the 
other case |T(ed)|. This prescription along with the nor- 
malization condition fixes the weight of the continuation 
process. 

In vertices C, D and G every weight for the exit path 
hopping to or out of the impurity is chosen as 1x14/2 
so that the forward and backward probabilities in G are 
equal. Continuation and bounce weights in vertices C 
and D are determined in the same way as for vertex B; 
the bounce weight turns out to be |rek| if forward contin- 
uation lowers the weights, otherwise zero. In G, there is 
no bounce back. In vertex F a hopping from site k to site 
q occurs with weig 

^ rain[\v^\,\v,\] ^ Finally, as discussed 
earlier, in vertex E the path is always forced to continue. 

The above rules can easily be extended to the case in 
which the directed loop enters a vertex on an empty site. 
Although the above rules satisfy detailed balance, they 
are in no way unique. We chose the above rules after 
testing a few other possibilities since they were similar 
in efficiency if not better than other ones. It has been 
suggested that if the bounce probability is large then the 
algorithm is likely to become inefficient, since the pro- 
posed change is being rejected. In our case a large bounce 
weight in electron sites far away from the Fermi surface 
is natural since the electron sites there either remain oc- 
cupied or empty most of the time. On the other hand, 
a large or U leads to a large bounce probability in 
the vertex B at the impurity. Thus, a large or U may 
cause inefficiencies in our algorithm. One can choose a 
different set of weights to reduce the bounce weights by 
taking into account e/c at all values of fc, but determining 
the weights in this case becomes difficult. We have tested 
a few different set of weights which gives less bounce back 
at the impurity and band electrons. Unfortunately, our 
attempts have not improved the efficiency of the algo- 
rithm further. So here we report on the results using 
simplest algorithm discussed above. 

Until now the limit of r — > was not taken. As r is 
taken to the zero, one gets the continuous time version. 
From Aqk, Adk, and the bounce weights, one can easily 



evaluate the decay rates for the continuous time simu- 
lation. The Monte Carlo simulation in continuous time 
proceeds as follows: First, we pick a starting time, spin, 
and path direction. For the starting site, the impurity 
site is picked more often. Typically 50% of the starting 
points are at the impurity and the remaining 50% on the 
levels in the band with equal probability llq. The path 
for the loop continues in time until a decay occurs into 
one of the possible vertices. Then, the new level and the 
direction are determined by the exit process. If a path 
hits a time slice where the configuration changes before 
a decay occurs, then the vertex at that time slice is used 
to decide the exit process. The loop update continues 
until it closes. As the loop is constructed the occupation 
states along the loop are flipped. 



III. THE MULTILEVEL ALGORITHM 

In a given configuration C electrons hop between the 
band and the impurity site so that in a periodic config- 
uration in imaginary time, the electrons permute their 
positions. Due to the Pauli principle, this causes config- 
urations to have a positive or a negative sign. This is 
the reason for the factor a[C] in Eq.Q. Any physical 
quantity O can be computed using 



(O) = -TrOe 



Ec Oa[C]W{C) _ (Oa) 



Ec <^[C] W{C) 



(9) 



where the final expectation values are computed using 
the Monte Carlo algorithm discussed in the previous 
section that generates configurations with probability 
P{C) defined in Eq.(|SJ). Unfortunately, as the temper- 
ature decreases, both the numerator and the denomina- 
tor decrease exponentially which makes the calculations 
of fermionic observables at low temperatures extremely 
difficult. Thus one needs an efficient method to compute 
exponentially small numbers by averaging large positive 
and negative numbers, a problem that is generically re- 
ferred to as the Sign Problem. 

Recently, a clever trick referred to as the multi-level 
algorithm was discovered in the context of lattice QCD to 
compute exponentially small numbers [T9l | . In particular 
it was possible to compute the potential V{R) between 
quarks and anti-quarks, by computing the corresponding 
exponentially small Boltzmann weight exp(— y(i?)/T) at 
a temperature T. In lattice QCD this quantity can be 
computed by averaging the Wilson loop which is typically 
of order 1 for a given configuration. It was shown that the 
multi-level algorithm could compute averages of Wilson 
loops that were as small as 10~^°. The basic idea was to 
write the observable, in this case the Wilson loop, as a 
product of many terms such that each of the terms is not 
very small even though the product is very small. In this 
article we show that a similar approach can be applied to 
compute the average sign in the fermionic problem using 
the directed-loop algorithm. 



5 



In order to apply the multi-level algorithm let us divide 
the Euclidean time j3 of the lattice into 2^^^ parts with 
the same time thickness. Let us denote the sub-lattice 
configuration of fermions inside each of these parts by 
Csiji = 1,2, ...,2^~^ and the boundaries between the 
sub-lattices by Cb^ where r = 0,e,2e, f3 represents 
the Euclidean times at the boundaries. Periodic bound- 



ary conditions in time means that C;,, 



C, 



Now, 



the boundary configurations Cf, and sub-lattice config- 
urations Cs determine the entire configuration C. The 
probability P{C) can then be expressed as 



(10) 



i=l 



where P(Cb) is the probability of finding the configura- 
tion Cb on the boundaries and P{Cb,Csi) is the condi- 
tional probability of finding the configuration Cg. given 
the boundary configurations Cb- Clearly, P{Cb,Csi) de- 
pends only on the boundaries that bound the ith sub- 
lattice. Since the sign of a configuration a[C] can be 
written as a product of sign factors coming from each of 
the sub-lattices, the average sign can be written as 

{a)= P{Cb)l[a,P{Cb,Cs^), (11) 

where at is the sign that comes from the sub-lattice i. 
Now J2c '^i P{C^b,Csi) is just the average sign of the 
i-th sub-lattice with a fixed boundary configuration. So 



{a) =J2 P{Cb)Y[{a,){Cb), 



(12) 



where {ai){Cb) is the average sign of the sub-lattice i with 
the boundary Cb- 

The multi-level algorithm proceeds as follows: First, 
to generate a sequence of boundaries Cb, one updates 
the entire lattice. Then, with a fixed boundary config- 
uration Cb, one generates a subsequence of Ns config- 
urations for each sub-lattice. The directed loop algo- 
rithm is well suited for this update; when the directed- 
loop encounters the fixed time slice the path is forced 
to bounce back! Clearly one can estimate the average 
sign ((Ti) of each sub-lattice independently using the Ng 
configurations and then use their product to compute 
(cr). It should be emphasized that although the sub- 
lattice averages are not exact there is no systemic errors 
in this approach. Ns is determined empirically so as to 
make the calculation efficient. The sub-averages do not 
have to be calculated more accurately than the size of 
the fluctuations due to change in the boundaries. Once 
the sub-lattices are updated the entire configuration is 
again updated to generate a new set of boundary con- 
figurations Cb- Repeating this process, a series of sign 
measurements are generated. The final result of the sign 
is obtained by averaging these measurements. The sta- 
tistical noise in the sign is reduced because effectively one 



FIG. 3: A schematic description of the four level algorithm. 



is summing over 
measurements. 



configurations for each of these 



In the multi-level algorithm one can also build nested 
levels of update. In particular one can performs a K 
level update if the Euclidean time direction is divided 
into 2^~^ parts. During the update at the first level 
all time slices t = Q,e, ...(2^~^ — l)e are held fixed, in 
the update as discussed above. At this level Ns loop 
updates are performed between the fixed time-slices. At 
the second level only time slices &t t = 0, 2e,4e,... are 
held fixed. After every update at the second level, the 
full first level update is performed. This procedure is 
repeated Ns times. This process is repeated at higher 
levels; such that at each highest level update, all the 
complete lower level updates are performed. Thus, at 
the Kt\i level the whole lattice is updated with just the 
t = time-slice held fixed. One can extend the ideas of 
measuring the sign for a single level update as discussed 
above, to the K-\eve\ update. After summing over the 
signs of all the configurations generated during the K 
level update one computes a value of the sign. If at each 
level one performs Ns updates for each sub-lattice, one 
performs on the order of {Ns)^ x 2^~^ loop updates after 
the Kth-\eve\ update. Since /? = 2^~^e, then for a fixed 
e if the effort for a single loop update remains fixed, the 
effort to compute the sign for a full K-\cvc\ update grows 
as a power of /?. As more levels are introduced, it takes 
longer for the measurement of sign. In Fig. 13 we show a 
schematic description of the multi-level idea. 

One can apply the multi-level technique to other ob- 
servables which can be written as product of quantities 
on each of the sub-lattices. We call such observables as 
being compatible with the multi-level algorithm. The 
optimum number of levels should be determined empiri- 
cally. This number can depend on the observable to be 
calculated. We also found that the full multi-level algo- 
rithm is the most efficient for fcrmion sign problems due 
to the large oscillations. 
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IV. OBSERVABLES 

Since the Monte Carlo update is performed in the occu- 
pation number basis, aU diagonal observables 0[n] that 
are functions of the occupation numbers can easily be 
calculated using the formula 



(0) = 



(13) 



Average occupation number of a level is one of the ob- 
servables which belongs to this class. Another important 
diagonal observable in the Anderson impurity model is 
the local susceptibility 



Xlo 



(14) 



which can be obtained from configurations of the impu- 
rity as a function of imaginary time. This method of 
computation can only be reliable for regions of tempera- 
ture where the sign problem is mild or when the multi- 
level algorithm discussed above is applicable and useful. 
Both, the average occupation numbers and the local sus- 
ceptibility, are compatible with the multi-level algorithm 
and hence the algorithm can be used to alleviate the sign 
problem in calculating these quantities. 

The quantity more directly relevant to experiment is 
the impurity susceptibility, Xim, which is the total sus- 
ceptibility minus the free susceptibility: 



Xtot 



Xi 



where Xi is the susceptibility from the i-th free electron 
site. In the determinantal Monte Carlo (the Hirsch and 
Fye algorithm) the impurity susceptibility is very difficult 
to calculate due to statistical noises in the simulation 
The problem is that one has to calculate the to- 
tal susceptibility for the Anderson Hamiltonian and then 
subtract the susceptibility for the free case from it. The 
total susceptibility is a quantity of order N (the number 
of band electron sites) , but the impurity susceptibility is 
of order 1. So one has to calculate the total susceptibil- 
ity with the error of order 1 or less, which is extremely 
difficult for large lattice N . From the Clogston- Anderson 
compensation theorem |ll) . for a large bandwidth with 
a flat energy density and equal hopping amplitudes, one 
can expect that the local susceptibility is equal to the 
impurity susceptibility. But in the study of mesoscopic 
fluctuations the two susceptibilities can be quite differ- 
ent [221 . One of the main advantages of our method is 
that with our algorithm the impurity susceptibility can 
be measured with significantly reduced statistical noise. 
Let us now discuss how we can compute Xim- 

For configurations in the occupation number basis gen- 
erated during the update, the hopping of electrons oc- 
curs at only a small number of electron sites. The rest 
of the energy sites appear to be free; an advantage of 



working in the "momentum" space lattice. Suppose that 
Miop sites are involved in the hopping for a given con- 
figuration C . Denote the configuration of those sites by 
Chop, and the free part by Cf. Then, one can express the 
probability of C as P(C) = P(Chop)P(Chop, C/), where 
-P(Chop,C/) is the probability of Cf with a given Chop. 
The impurity susceptibility can be expressed with P{Cf) 
and P(Chop, Cf) as 

Xim X! ^(Chop)CT[Chop] 

C'hop 

X {x(Chop) 

+ E^(^hop,C/)[x(C/)- ^ X.- E ^']}' 



Cf 



where x(C'hop), x{Cf) are the susceptibilities from sites 
that contain electron hops and those that appear free in 
a given configuration. Since C/ contains no hop, one can 
see that 

5]p(Chop,C/)[x(C/)- 5]x.] = o. 

Using this one can find that the impurity susceptibility 
is given by 

Xim = -7^ /cr[Chop][x(Chop) - ^*]/' 

where the free susceptibility of the sites in Chop only 
needs to be subtracted. In this method the size of sta- 
tistical fiuctuation of the measurements of the impurity 
susceptibility is of order of the number of electron sites 
in which the hopping occurs for a configuration. So the 
statistical noise in this method are substantially reduced. 
The observables that go into Ea. (|15l) are again compati- 
ble with the multi-level algorithm. 

Unfortunately, we have found that the above technique 
is still noisy in practice. Interestingly, one can reduce the 
statistical noise in x(C'hop) further by using the technique 
of improved estimator that is commonly used in cluster 
algorithms. In this method one identifies all spin clus- 
ters for a given configuration Chop that can be flipped 
independently and performs a partial average over these 
cluster flips. Unfortunately, this step is not compatible 
with the multi-level algorithm. But for the case where 
the sign problem is moderate the impurity susceptibility 
can be computed very efficiently using our algorithm. 



V. RESULTS 

In this section we discuss results from the simulation 
of the Anderson Hamiltonian given in Eq.0J with the 
algorithm described in the previous sections. First, let 
us focus on the calculation of the average sign using 
the multi-level algorithm. For this purpose we choose 
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FIG. 4: Average signs versus /3 for iV = 750, D = 5, U = 2, 
Ed = -1 and r = 0.5. 



= 750 equally spaced energy levels with a band- 
width of 21) = 10. We choose Vk = V such that 
r = npV^ = 0.5 and U — 2. We have studied three dif- 
ferent temperatures by choosing /3 to be 40, 80 and 160 in 
order to see the effectiveness of the multi-level algorithm 
in the computation of the average sign. The sub-lattice 
thickness e is chosen to be 10 so that at /3 = 40 we have 
four, at /3 = 80 we have eight and at P — 160 contains 
sixteen sub-lattices. 

For the sub-lattice of thickness e = 10, we found that 
Ns = 10 updates was necessary to get an reasonable esti- 
mate of the average sign of the sub-lattice. To complete 
one full cycle of the all multi-level updates, (2iVs)^/2 
loop updates are required where K = 3 for f3 = AO, K ~ 4 
at /3 = 80 and K — b a± (3 — 160. In addition, at each 
higher level sub-lattices were updated 4 times to generate 
new boundaries between sub-lattices. 

At /3 = 40 and 80, we computed (cr) to be 4.99(11) x 
10~^ and 4.13(16) x lO"'', respectively, where the errors 
are of the order of a few percent. At /3 = 160 the sign 
average is so small and the projected time is so long that 
the simulation was stopped when the error was about 30 
percent. The average sign at /3 = 160 was 4.7(1.2) x 10^®. 
The computational time taken for these results are 3hrs., 
92hrs., and 2000hrs., respectively. For the /3 = 160, 6 
CPUs were used with different random number seeds to 
collect the data. In Fig. 01 we plot the average sign as a 
function of (3 and we see that all the values fall nicely on 
an exponential form as expected. 

Now the biggest question to answer is whether the 
multi-level algorithm is useful. We would first like to 
point out that it is still difficult to compute the average 
sign with reasonable errors at very small temperatures. 
The required effort grows as at least a large power of 
/3, and we cannot rule out an exponential growth at the 
moment. However, the fact that for /3 — 160 we could 
compute numbers of the order of 10~^ itself is an indi- 
cation that some progress has been achieved. Without 
the multi-level algorithm this would have been impossi- 
ble. If we look at the individual values of the sign com- 
puted by the multi-level algorithm after each update we 
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FIG. 5: Monte Carlo sequence of fermion signs for /3 = 40, 
and 160 from the top. 
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FIG. 6; The first 2000 of the positive(left) and negative signs 
in logarithmic scale of the (3 — 160 run. 



learn something further. Fig. |S1 shows these values of the 
signs for different values of f3. We see that using the 
multi-level algorithm the values of the signs are domi- 
nated by positive values. For example one can see that 
at /3 = 160 the average sign has most contributions from 
the positive side. We show this in Fig. by focusing on 
the first 2000 positive and negative values. We see that 
although the very small values come with equal weight 
between positive and negative values, there are a few pos- 
itive numbers that dominate over the negative numbers. 
The multi-level algorithm has allowed us to find config- 
urations whose signs, when summed up, leads to a very 
small number which can either be positive or negative but 
occasionally it also leads to numbers which are orders of 
magnitude larger but mostly positive. The large error in 
the average sign is mainly due to these large fluctuations 
but always in the positive direction. In the next section 
we will discuss why this observation is interesting. 

In order to show that the new algorithm is indeed in- 
teresting, for the parameters chosen above we compute 
the average occupation numbers for the various energy 
levels. In Fig. [7|we plot the differences between the av- 
erage occupation numbers in the interacting cases and 
the free cases at the temperatures /3 = 25 and 50. The 
open circles give the results in the cases where ?7 = in 
which case one can obtain the result also from the exact 
Green functions assuming the density of energy levels is 
smooth (solid line). Thus, we see that our algorithm can 
indeed produce results in good agreement. For the case of 
U = 2, the Kondo temperature for this system is roughly 
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FIG. 7: Occupation numbers minus the free Fermi function 
at /3 = 25(top) and 50(bottom) for iV = 750, 73 = 5, ed = -1, 
and r = 0.5. 



Tk = 0.09 as seen from the numerical renormalization 
method 0, 0| ■ In this case (3 — 25 and 50 correspond 
to a temperature of about Tfc/2 and Tk/4:. We see that 
indeed the Kondo resonance appears as expected. 

Finally, we focus on the local and impurity suscepti- 
bilities. To check the Clogston- Anderson compensation 
theorem, a larger bandwidth 2D = 20 is chosen with 
TV = 2000. For [/ = 2(ed = -1), Xio and Xim are shown 
in Fig. IHl For the local susceptibility, with the multi- 
level method we were able to calculate the local suscep- 
tibility at a much lower temperature than the impurity 
susceptibility. We find that xio and Xim are in reason- 
able agreement as expected from the Clogston- Anderson 
compensation theorem. We also compare our result with 
the NRG curve (solid line) obtained for Tk ~ 0.08 which 




FIG. 8: Local and impurity susceptibilities for U — 2(top) 
and U = 25(bottom) with U/T = 4 fixed. 



passes through most of the data points. 

Since our simulations are in continuous Euclidean time, 
we can simulate a large U without increasing the dis- 
cretization error. In the limit of large U, the Ander- 
son model converges to the Kondo model. In the Kondo 
model, band electrons and impurity spin interact with 
the coupling J. From the Schrieffer- Wolff transformation 
[23], the effective coupling J of the Anderson Haniilto- 
nian for a large U is In order to go towards the 

Kondo limit wc fix U/T — 4 and study the case where 
U — 25{€d — —12.5). The local and impurity suscepti- 
bilities are plotted in Fig. |S1 We see that now these two 
are completely different. We attribute this difference to 
the fact that U ^ D in which case the compensation 
theorem is no longer valid. 

We have also computed the various observables dis- 
cussed in this article for the Hamiltonian that contains 
mesoscopic fluctuations. We find that typically we can 
use our new method to compute quantities for tempera- 
tures as low as r ~ 7fe/4. Using the multi-level technique 
we can also go down to temperatures of about T ^ Tk/10. 



VI. CONCLUSION 

In this article we have investigated a new algorithm for 
a model involving a band of fermions interacting with a 
single impurity in the occupation number basis in "mo- 
mentum" space. We use the efficient directed loop algo- 
rithm to update the configurations and absorb the sign 
into observables. We find that the sign problem is mild 
down to temperatures of order Tfe/4. Further, the new 
approach allows us to explore a new multi-level algorithm 
to compute average signs efficiently. We were able to 
compute signs of the order of 10~* with moderate ef- 
fort. The new approach also allows us to calculate certain 
quantities like the impurity susceptibility more efficiently 
than conventional Monte Carlo methods. Finally, our al- 
gorithm can easily be extended to several impurities. 

The average of the sign over configurations that are 
generated in the multi-level algorithm fluctuates between 
small values which can be both positive and negative and 
large values which are orders of magnitude larger but al- 
ways positive. The effort for this grows as {2Ns)^ /2 
where f3 = 2^~^. Although this does not solve the sign 
problem completely, since the positive numbers can still 
fluctuate a lot, perhaps half of the sign problem has been 
solved. If this is true then we think this is an excit- 
ing step in the solution to the full sign problem based on 
the recent progress in solving certain sign problems using 
the meron cluster algorithm Q . There it was possible to 
rewrite the partition function in terms of configurations 
where the Boltzmann weight was either zero or positive. 
Thus all negative signs were eliminated. The second step 
was algorithmic when all zero configurations were elimi- 
nated in an accept reject step. An intriguing question is 
whether something similar can be achieved in the present 
case. We leave this question for future research. 
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